Why is Random Close Packing Reproducible? 
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We link the thermodynamics of colloidal suspensions to the statistics of regular and random packings. Ran- 
dom close packing has defied a rigorous definition yet, in three dimensions, there is near universal agreement 



on the volume fraction at which it occurs. We conjecture that the common value of < 



0.64 arises from a 



divergence in the rate at which accessible states disappear. We relate this rate to the equation of state for a hard 
sphere fluid on a metastable, non-crystalline branch. 



When Bernal turned to sphere packings as a route towards 
understanding liquids [ II], he recognized the importance of the 
random close-packing density, <p Icp . This density was defined 
operationally to be the fraction of occupied volume in the 
densest disordered packing of hard spheres. Early measure- 
ments of (f) Icp by Scott [2] and Bernal and Mason [3] involved 
pouring, shaking and kneading ball bearings in flasks and bal- 
loons. Since then, these measurements have been reproduced 
by countless experiments and numerical algorithms, which 
find a common value of </> rcp w 0.64 in three dimensions fl. 
There are therefore many distinct but numerically-consistent 
definitions of rcp that depend on procedure. At the time 
of Scott's experiments, Bernal and Mason wrote, "The fig- 
ure for the occupied volume of random close packing-0.64- 
must be mathematically determinable, although so far as we 
know undetermined." [3]. Nearly fifty years later, however, a 
mathematically-rigorous definition of random close-packing 
still remains elusive. Why should an ill-defined state have a 
reproducibly accepted value for its volume fraction? 

Two approaches towards a mathematically-rigorous defi- 
nition have been suggested recently. Torquato, et al. ^ 
pointed out that there is an inherent discrepancy in the con- 
cept of a "random close-packed state," because one can al- 
ways obtain a more closely-packed state (a denser state) by 
introducing order into the system. This must be true be- 
cause 4> lcp is lower than the maximum close-packed den- 
sity Bg], </>fcc = VT/3V2 ~ 0.74, corresponding to the close- 
packed FCC crystal. Torquato, et al. [51] proposed an alter- 
nate way of thinking about random close-packing, in terms of 
a "maximally-random jammed" (MRJ) state at 4> rcp . "Maxi- 
mally random" states or configurations are those with minimal 
values of typical order parameters, such as bond-orientational 
order or crystalline order. "Jammed" states have the property 
that any particle or set of particles cannot be translated with 
respect to any of the rest of the particles in the system with- 
out introducing overlaps. Kansal, et aJ. |7[] showed that sev- 
eral different order parameters yielded consistent estimates for 
the packing fraction of the maximally random jammed state, 
0mrj ~ 0.63, in good agreement with </> rcp . This is very 
encouraging. Nonetheless, there is some uncertainty in this 
approach associated with the order parameter; one cannot cal- 
culate all possible order parameters, or determine that all order 
parameters yield the same MRJ packing fraction. Moreover, 
there may be a specific type of order associated with the MRJ 
state lllll . which would then be maximal rather than 

minimal in the MRJ state, though these special metrics are 



likely to be of a different nature from those studied in ^ . 

A second approach by O'Hern, et al. is based on the energy 
landscape of systems of soft spheres. Specifically, Refs. Bill 
[l2ll considered potentials of the form 



V(r) 



r/cr) a /a 




for r < a 
for r > a 



(1) 



Here, e is the characteristic energy of interaction and a is the 
particle diameter. For a — 3/2, a — 2 and a = 5/2, O'Hern, 
et al. studied the fraction of ideal gas states fo(4>) that belong 
to basins of attraction corresponding to zero energy states (or 
equivalently, allowed hard sphere states). It was found that 
— dfo((j>) / dip, i.e. the rate at which the fraction of ideal gas 
states belonging to basins of attraction of hard sphere states 
shrinks with increasing density, appears to develop a delta- 
function peak at (f> rcp in the infinite system size limit lfTll[T2Tl . 
Thus, they suggest that in the thermodynamic limit, the vast 
majority of ideal gas states belong to basins of attraction of 
hard sphere states that jam at <j> lcp . Here, a state is "jammed" 
if there are no zero frequency vibrational modes except for 
those due to floaters (particles with no overlapping neighbors) 
or overall translation and rotation of the system 11311 . The 
uncertainty in this approach is associated with how the energy 
landscape is explored 11141 fl5ll ; different algorithms may yield 
different final energy minima for a given ideal gas state. 

In this paper, we exploit a relation that has been used to 
calculate the pressure of a system to explore the connection 
between the pressure - a thermodynamic quantity that can be 
measured for colloidal suspensions - and the behavior of hard 
sphere packings. In particular, we show that the fractional 
rate at which allowed states disappear with increasing volume 
fraction is proportional to the pressure. We also show that free 
volume theory provides a reasonable fit to the equation of state 
of the hard sphere liquid, with a single fit parameter 4> max , 
corresponding to the density at which the pressure diverges. 
We find that max = 0.640±. 006, in excellent agreement with 
4>rcp- We conclude by discussing these results in the context 
of previous work and making some conjectures regarding the 
origin of the reproducibility of </> rcp . 

We first derive a useful relation between the pressure and 
the rate at which allowed states disappear. This relation has 
been used to calculate the pressure from a Monte Carlo run 
for hard particles by calculating the minimum density change 
needed to introduce overlaps between particles 11611 . Here, 
we use the relation to gain insight into the packings of hard 
spheres. We consider the number of available states as a func- 
tion of volume fraction <f> — pv, where p is the number density 
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of spheres and v is their volume. We consider the probability 
of finding an allowed hard sphere state at p, 



K(p) = Z HS (p)/Z 1G (p), 



(2) 



where Zhs and Z\q are the numbers of allowed configura- 
tions (or equivalently, the canonical partition functions) for 
hard spheres and for an ideal gas of point particles, respec- 
tively. We now consider the effect of increasing the vol- 
ume fraction of a particular state at density p. The proba- 
bility of not finding a new state at p + Sp, given this state 
at p, is J(p)6p, where J(p) — Using the thermo- 

dynamic definition of pressure, p = Tdln Z/dV, we find 
J (p) = W [phs^) - Pig(p)] > where p m and p IG are the 
pressures of the hard sphere and ideal gas systems, respec- 
tively, V is the volume and T is the temperature (the Boltz- 
mann constant is unity). Equivalently, 



J{4>) = [phs W - 4T/v] . 



(3) 



In Fig. 1(a) we show the experimentally-measured equa- 
tion of state for a hard sphere system (circles) fl7tl . The 
pressure increases with <f> until <px w 0.49, at which point 
fluid and crystal begin to coexist. At <f> w 0.54 the system 
is completely crystalline and phs(</>) eventually diverges at 
FCC = 7I-/3V2 ~ 0.74, the packing fraction of the close- 
packed FCC lattice. Eq. (O implies that J(<f) also diverges 
at </>fcc> signaling an inability to construct any higher density 
states, i.e. close packing. Thus, if all the states are counted in 
Zhs, J does not diverge until 4>fcc- 

There is also a metastable branch of the equation of state 
that is reproducibly measurable numerically (see Fig. 1) 
lfl8l[l9ll and potentially measurable experimentally in rapidly 
sedimented colloidal suspensions ll2b1l . This branch is used to 
study the colloidal glass transition 112 ill and dynamics in gran- 
ular me dia 12211 . Along this branch, the pressure apparently 
divergesil8l ll9|] at rcp as pns(4>) ~ (</>rc P - Such a 

divergence is predicted on the basis of polytope theory, or free 
volume theory, for classical hard spheres |H,|23t]. According 
to Eq. (01, this would lead to a divergence in J there. To 
construct an analytical approximation to a metastable branch 
that continues smoothly from </>x to a divergence at some 
0max, we therefore turn to free volume theory. To calculate 
the free energy within free volume theory, we must construct 
the Voronoi tesselation for any particular packing ll24ll . For 
a given packing, each sphere is allowed to independently ex- 
plore the free volume of its cell. The free energy involves a 
sum over all possible tesselations and requires knowledge of 
the distribution of free volumes lf25ll . However, the change 
of free energy with <fi, or the pressure, should not depend on 
the tesselation sufficiently close to </> max . If 4> is decreased in- 
finitesimally near max , the increase in free volume of each 
cell for an isotropic, affine rarefaction of the packing scales as 

[(^max/^) 1 ^ 3 — l] as long as the shape of each free volume 
region remains fixed. Thus, near </> max , the pressure depends 
on only one free parameter, <p max : 
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FIG. 1 : (Color Online) Measured equation of state p (in units of pT) 
versus volume fraction <f> for a hard sphere fluid from experiments by 
Rutgers, et al. 11711 (circles) and simulations by Rintoul and Torquato 
1 1 811 (diamonds) and Speedy [19] (squares). The dashed curve is the 
Carnahan-Starling equation of state. In a) we show two fits to free 
volume theory, corresponding to m ax = 0.635 (solid), which is the 
best fit to the experimental data for the equilibrium liquid branch, 
and <^ max = 0.6465 (dotted), which is the best fit to the numerical 
data for the metastable branch. In b) we show free volume theory 
for <f> max = 0.6465 (dotted) and max = 0.640 (solid). Part (a) of 
the figure shows the quality of the fit to the equilibrium liquid branch 
while part (b) shows the quality of the fit to the metastable branch. 
Inset to (a): the variation in the mean-square error A between free 
volume theory and the experimental data as a function of (/> m a X - 



PFV 



T<t> 2 d 



(4) 



where D is the dimension of space |24]. Note that as the vol- 
ume fraction further decreases, the geometry of the allowed 
volume will change and, though the scaling factor remains the 
same, the prefactor will change, leading to additional correc- 
tions. Recent arguments [26] based on statistical geometry 
suggest that Eq. should be at least approximately correct 
for the metastable branch. 

Fig. 1(a) shows that free volume theory provides a very 
reasonable fit to the equilibrium liquid branch of the data of 
Ref. lfl7ll just below fix, with (f) max ~ 0.64 in 3D, in excel- 
lent agreement with measured RCP values. The inset shows 
that A, the mean-squared error, is a strong function of (/> max . 
We find a comparable value of max w 0.636 when we fit 
the range <fi £ [0.3, 0.5] to the first ten virial coefficients |f27 
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Unsurprisingly, we also find a similar value (</> max ~ 0.635) 
when we fit to the Carnahan-Starling approximate equation of 
state lElll . The latter function is known to describe experi- 
mental or numerical measurements of the liquid branch of the 
equation of state to within measurement error lfl7[|29i[30l[3Tll . 
It is worth noting that the fitted values of max are insensi- 
tive to the region over which we fit the data. Using the entire 
equilibrium liquid branch (fi £ [0, 0.5] only changes </> max in 
the last decimal place when fitting to the experimental data, 
the virial expansion, or the Carnahan-Starling approximation 
(for example, we obtain max « 0.637 when fitting to the 
Carnahan-Starling equation over the whole range). Aste and 
Coniglio B32I1 have used local random configurations of hard 
spheres to estimate the pressure in free volume theory and 
there, too, find an independent branch of the equation of state 
diverging at max m 0.65. 

We have also fit free volume theory to numerical measure- 
ments of the metastable branch Il8lll9ll . Here we find that the 
best fit corresponds to <f> max ~ 0.6465-less than 2% higher 
than the result for the fit to the equilibrium liquid branch. In 
Fig. 1(a), we plot Eq. (4) for both max = 0.635 (solid) and 
4>max — 0.646 (dotted); these curves are nearly indistinguish- 
able in terms of their ability to fit the equilibrium portion of 
the curve. In Fig. 1(b), we plot the same data and fits over a 
narrower range of (fi and larger range of pressure to show the 
quality of the fit to the metastable branch. Note that our one- 
parameter fit reproduces not only the shape of the divergence 
but the overall amplitude of the pressure reasonably well. The 
final value <fi max = 0.640 ± 0.006 represents the range of fit 
values that we obtain, though the metastable branch strongly 
favors the high end of this range, as seen in Fig. 1(b). 

We now note that the many different algorithms and prepa- 
ration methods that yield rcp « 0.64, have one point in 
common: they all contrive to avoid crystallization and should 
thus correspond to a metastable branch of the pressure. For 
two very different approaches, conjugate gradient minimiza- 
tion llllLll2h and the Lubachevsky-Stillinger algorithm SB, 
there is good evidence that the accessible states follow the 
branch estimated in Fig. 1. Both procedures yield a pair cor- 
relation function g(r) that diverges at contact as \(fi rC p — </>| -1 , 
implying a likewise diverging pressure. These approaches 
both observe a catastrophic loss of states, i.e. they jam, at 
(fi Tcp , consistent with Eq. (1). 

These results suggest that (fi rcp corresponds to a divergent 
endpoint of a metastable branch of the equation of state for 
hard spheres. This provides another way of thinking about 
why random close-packing has not been a well-defined con- 
cept It is well understood that metastable branches are 
somewhat arbitrary. To obtain metastable branches, one must 
introduce additional constraints that effectively exclude cer- 
tain states (such as crystalline states) from the partition func- 
tion rf34ll . Different choices of constraints lead to different 
metastable branches. Alternatively, one can think in terms of 
the Andreev-Fisher droplet model l35i l36ll . or the instanton 
approach to first-order transitions [37], which predict an es- 
sential singularity at the onset of crystallization, (fix- This es- 
sential singularity precludes analytic continuation of the pres- 
sure beyond (fix', physically, droplets of the nucleating crys- 



talline phase prevent the clean definition of the metastable 
branch. 

The possibility of multiple metastable branches leads to the 
possibility of multiple divergent endpoints. The existence of 
multiple divergent endpoints is strongly suggested by results 
of Torquato, et ai. la*] using the Lubachevsky-Stillinger algo- 
rithm [33]. In this algorithm, one starts with an equilibrium 
liquid configuration at low volume fraction and compresses 
by increasing the diameters of all the particles at some rate T. 
The system jams at some packing fraction <fif(T), which ap- 
proaches 4> rcp from above in the limit T — > 00. This suggests 
that one can find metastable branches that end at any value of 
(fi between tfi lcp and <^>fcc ~ 0.74. 

Torquato, et ai. have suggested that the state at (fi Tcp is 
special in the sense that it corresponds to a "maximally ran- 
dom jammed" state la]. It seems reasonable that this con- 
cept should be generalizable to an entire "maximally-random 
jammed" metastable branch of the pressure that ends at </> rcp . 
This is the philosophy underlying calculations of Rintoul and 
Torquato IU8I1 . who followed a metastable branch of the pres- 
sure by discarding all states with appreciable values of the 
bond orientational order parameter. In that case, they found a 
diverging pressure at (fi TC p- 

Another possibility that is not inconsistent with the scenario 
of Torquato, et aJ. is that (fi lcp might be better defined than 
the metastable branch that it terminates. It might correspond 
to a true singularity of the free energy that is inaccessible in 
equilibrium due to the essential singularity at (fix, the onset 
of crystallization. Exact analyses of one-dimensional mod- 
els show that it is possible for a system to have a metastable 
branch that is not well-defined but that ends in a divergent 
endpoint that is a true singularity B38I1 . 

We therefore conjecture that (fi lcp represents a special well- 
defined divergent endpoint of a set of metastable branches of 
the pressure. Any procedure that samples a nonzero fraction 
of states belonging to metastable branches with this endpoint 
will yield a divergent pressure, and therefore a divergent rate 
at which states disappear, at rcp . If this conjecture is cor- 
rect, it would explain why so many different procedures, all 
sampling states somewhat differently, yield the same value of 
(fi Tcp . At the same time, it is clear that other procedures might 
yield different jamming densities by avoiding states that be- 
long to metastable branches with endpoints at (fi Tcp . 

Based on results for hard spheres obtained from replica the- 
ory, Parisi and Zamponi [39, 40] have suggested a similar but 
more elaborate scenario. In this picture, (fi rcp is the divergent 
endpoint of a metastable branch, but there is also an another 
point on the branch, <fi g < <fi rcp , which marks a thermody- 
namic glass transition. Above (fi g , the configurational entropy 
vanishes; the system must remain in the lowest free energy 
states and the relaxation time is infinite. This scenario is not 
inconsistent with our conjecture. 

A singularity in the rate of change of number of states does 
not necessarily imply that most initial states will flow to <^> rcp . 
It is possible that most states might have their jamming thresh- 
olds at values of (fi below <fi lcp , leaving only a few that termi- 
nate at (fi lcp . Numerical results suggest that the opposite might 
be true. Indeed, they suggest that (fi lcp may not only mark a 
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well-defined divergent endpoint of a set of metastable equa- 
tions of state, but that an even stronger condition might hold: 
the distribution of jamming thresholds may actually have a di- 
vergent maximum at <fi lcp in the thermodynamic limit. This 
conjecture is motivated by results of O'Hern, et al. jTll,[l2tl . 
which suggest that for several soft repulsive potentials, the 
overwhelming number of ideal gas states belong to basins of 
attraction of hard sphere states that have their jamming thresh- 
olds at 4> lcp . Note that the probability of belonging to a basin 
of attraction of a state with a jamming threshold at tj> depends 
on both the size distribution of basins of attraction and the 
distribution of jamming thresholds. For small systems, Xu, et 
al. 114 111 have separated the these two distributions by explic- 
itly enumerating the jamming thresholds; their results suggest 
that the distribution of jamming thresholds is maximal near 
4> lcp , consistent with this conjecture. While it is unlikely that 
the maxima of the distributions of jamming thresholds and of 
jamming thresholds weighted by their basin of attractions and 
the divergent endpoint of the metastable branch all coincide 
at exactly the same density, it is possible that they agree with 
4> = 0.64 to within 1-2%, so that any one of these could con- 
stitute a reasonable definition of 4> rcp l42ll . 



The quantity lZ((f>) defined in Eq. (0, the ratio of hard 
sphere states to the total number of ideal gas states at a 
given <fi, includes not only hard sphere states at their jamming 
thresholds at <fi but also hard sphere states below their jam- 
ming thresholds. In this paper, we have argued that the form 
of the equilibrium equation of state suggests that at (f> lcp , the 
latter states completely dwarf the number of hard sphere states 
that are at their jamming thresholds, hiding the divergence in 
the pressure. It is only when these more ordered states are 
excluded by restricting the system to a metastable branch of 
the pressure that a signature of <fr rcp appears. This reasoning 
suggests that a clean definition of <fi lcp will rely not only on 
the distribution of jamming thresholds, but the number of al- 
lowed hard sphere states at their jamming thresholds at 4> rcp 
as a function of system size 04311 . 
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